Critical behavior of the three-dimensional compressible Ising antiferromagnet at 

constant volume: a Monte Carlo study 
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Extensive Monte Carlo simulations in the semi-grand-canonical ensemble are used to study the 
critical behavior of a three-dimensional compressible Ising model with antiferromagnetic interactions 
under constant volume conditions. Elastic forces between spins are introduced by the Stillinger- 
Weber potential and energy parameters are chosen in such a way that antiparallel spin ordering 
is favored, analogous to the antiferromagnetic coupling in the rigid Ising Hamiltonian. All the 
quantities analyzed strongly indicate that the system remains in the universality class of the standard 
(rigid) three-dimensional Ising model, in contrast with theoretical predictions. 
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I. INTRODUCTION 

The Ising model, first proposed by Lenz in 1925 as 
a microscopical model for ferromagnetism, has a (con- 
stant) interaction J included only between the nearest 
neighbors spins on a (rigid) lattice. Despite its (formal) 
simplicity an analytical solution has been found only for 
one-, and two-dimensional models in zero magnetic field i 
For the three-dimensional Ising model only approximate 
solutions are available (e.g. series expansions,^ renormal- 
ization group calculations^, e-expansionsf^ Monte Carlo 
renormalization group calculations^, and Monte Carlo 
(MC) simulations),^ although quite precise results have 
been found. 

It was clear, a few years after its introduction, that the 
Ising model could be employed to describe other phase 
transitions. For instance, phase separation in binary al- 
loys can be studied by an Ising model in which "up" and 
"down" spins are replaced with sites occupied by a "A" or 
"B"type atom, respectively. The chemical potential plays 
the role of the magnetic field, the density that of the 
magnetization, etc. and the appropriate statistical en- 
semble is the semi-grand-canonical, instead of the canon- 
ical. With the above mentioned analogy in mind, in the 
present work we will use the languages of magnetism and 
of alloys interchangeably. 

In a real magnetic crystal, however, atoms interact 
with each other through a combination of elastic and 
magnetic forces. The next step in realism is, therefore, 
the explicit introduction of elastic degrees of freedom in 
the traditionally rigid system. The resulting model is 
termed compressible Ising (CIM). Several empirical po- 
tentials have been proposed and employed to mimic the 
elastic force, from the simple Lennard- Jones to the more 
specific Tersoffj- Keating,^ Stillinger- Weber, ^ and oth- 
ers (these last three have been introduced, mainly to re- 
produce the interaction between silicon and germanium 
atoms in the study of Sii_a;Gea;, viewed as Ising binary 
alloy models). 

The issue of how the presence of elastic interactions 
affects the critical behavior of the Ising model has been 
intensively studied. In 1954, using thermodynamic con- 
siderations, RiceiS^ predicted that if an Ising system with 



divergent specific heat is put on a deformable lattice at 
constant pressure, it undergoes a first order transition. 
A number of calculations appeared lateriiitiSii^, all of 
which lead to the conclusion that a first order transition 
was expected to occur. All these studies assume that a 
"pure" Ising singularity happens at constant volume. In 
1968 Fisher— changed the situation drastically with his 
"hidden variable" theory, assuming that the "pure" phase 
transition occurs at fixed intensive variable, i.e. at fixed 
pressure. As a result of this hypothesis a second or- 
der transition was found with Fisher renormalized ex- 
ponents, if the critical exponent of the specific heat a 
is positive. A significant highlight in this controversy 
was the work of Larkin and Pikin^^ who, for the first 
time, considered a Hamiltonian with fluctuations of both 
the order parameter and the elastic modes and pointed 
out the special role of the macroscopic mode k=0 (in 
Fourier space). As a result a first order transition is 
found at constant pressure and quadratic coupling of the 
order parameter and the strain tensor, also if a = 0. 
It was later recognized that this result is only valid at 
low pressure, whereas at high pressure the critical be- 
havior seems to be much more complex— tii A different 
approach was used by Baker and Essaroi^ who mapped 
an Ising model on a compressible lattice, including the 
coupling between magnetic and elastic degrees of free- 
dom, onto a standard Ising (on a rigid lattice) , but with 
parameters which depend on the elastic degrees of free- 
dom. At constant volume the model exhibited identical 
critical behavior as the underlying rigid system, and not a 
first order transition as predicted. At constant pressure 
a second order transition was found with Fisher renor- 
malized exponents. This model, as well as others of the 
same type, however, considers a negligible shear modu- 
lus. This somewhat unphysical assumption was soon rec- 
ognized to be the reason of the disappearance of the first 
order transition. Further worlsiSiSSiSiiSS with more elabo- 
rate schemes did not change the result: a first order tran- 
sition was always predicted. Finally, we cite the long and 
complicated work by Bergman and Halperin^S^ in which 
a three-dimensional anisotropy in the elastic forces was 
introduced. In the isotropic case the same already estab- 
lished result was obtained, i.e. a first order transition at 
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constant pressure if a > 0, and a second order transition 
with renormalized exponents at constant volume. How- 
ever, this last condition is realized only if the atoms at the 
surface are fixed, otherwise the system has still enough 
degrees of freedom to develop a macroscopic instability. 
In the anisotropic case, however, a so-called microscopic 
instability was found at constant volume, and this is in- 
terpreted as a first order transition. In his Habilitation 
thesis Diinweg24 completely revised the topic. Starting 
from the Landau-Ginzburg- Wilson Hamiltonian he was 
able to show that the CIM under constant pressure ex- 
hibits mean-field critical behavior, in agreement with 
simulations At constant volume he predicted two 
first order lines ending in critical points, which are likely 
to belong to the mean- field universality class. If the mag- 
netic interactions are antiferromagnetic (AFM) instead 
of ferromagnetic's a quadratic coupling between the or- 
der parameter and the elastic deformation should be ex- 
pected. In this case the predictions are: (i) a second or- 
der phase transition with Fisher-renormalized exponents 
in the case of constant volume^S (this is in agreement 
with the e-expansion work of Bergman and HalperinSS) ; 
(ii) a first order phase transition in the case of antifer- 
romagnetic interactions and constant (zero) pressured 
As pressure increases, the first order line in pressure- 
temperature space should split into two first order lines 
at a triple point. The theoretical prediction in the case 
of FM interactions at constant volume has been checked 
by Tavazza et al.^^ by MC simulations. In disagreement 
with theory, they found a closed first order line which 
separates ordered and disordered phases. 

This result raises the intriguing question of whether or 
not theory will be correct for the case of AFM interac- 
tions and constant volume. This will be the topic of the 
study that we report here. 



II. MODEL AND SIMULATION TECHNIQUES 

We use the same model as considered previous in Ref. 
We consider a binary alloy of A and B atoms on 
the nodes of a distortable diamond lattice, free to move 
on the condition that the diamond four-fold coordina- 
tion is preserved, and that the atomic species on each 
node may change. The coordination requirement speeds 
up considerably the simulations as the list of the nearest 
neighbors of a given atom, which enter the Hamiltonian 
(see below), is known from the very beginning and does 
not change during the simulations. The diamond lattice 
is decomposed into eight interpenetrating simple cubic 
sublattices of linear size L, so there are N = 8L^ atoms 
(sites). Each atom in the system has four degrees of free- 
dom: the three spatial coordinate r and its species S, 
which is defined to be S" = 1 if the atom is an A, or 
S = —1 if it is a B type. The total number of atoms 
N is kept constant during the simulation, while the rel- 
ative concentrations of A and B particles can vary and 
are controlled by the chemical potential. The correspond- 



ing appropriate statistical ensemble is termed semi-grand 
canonical. The Hamiltonian is given by 

H = -^{HA - fJ.B)^Si + Hsw, (1) 

i 

where /ia, Mb are the chemical potentials of A and B, 
respectively, and Hsw is the Stillinger- Weber (SW) po- 
tential given by 

Hsw = H2bd + H^M- (2) 

The SW potential contains a two body interaction H2bd 
involving nearest neighbors and a three body i?3fcd inter- 
action that includes next-nearest neighbors as well. This 
last term is essential to stabilizing a diamond-like struc- 
ture. The two body interaction is given by 

H2M = S,)F2[r,j/Ro{S,,S,)], (3) 

where Rq is the ideal distance of the atoms, is the 
distance between sites i and j, F2 is a short range func- 
tion of the rescaled bond length (see Ref. |23for details), 
and e is the covalent binding energy. For a binary alloy 
of silicon and germanium (Si=A, Ge=B) it has been es- 
timated e(l, 1)=2.17 eV, e(-l,-l) = 1.93 eV, and e(l, -1) 
= 2.0427 eVrl and these values were used for of an elas- 
tic ferromagnet simulations^2i^ In the present work we 
increased the value for the A-B binding energy by 0.3 
eV to favor alternate ordering of A-B particles analogous 
to the antiferromagnetic ordering in magnets, specifically 
e(l, — 1) = 2.3427. The behavior for this value is expected 
to be typical of that in the AF regime, but a much larger 
value could conceivably produce unanticipated effects. A 
systematic study of the dependence upon the value of 
e(l, —1) is beyond the scope of this paper. The three 
body part of the Hamiltonian is given by 

HsM = {<S,,S,y/\{S,,Sk)'/^£{S,,S,,Sk) X 

F3[n,/RoiS,,S,),r,k/RoiSj,Sk)]icos9,jk + ^)'},(4) 

where F3 is a function of the same kind of i^2, is a 
simple function of the atomic species (see Ref. |23 for 
details), and cosSyfe = Tij ■ rjk/\rij ■ rjk\- The sum is 
performed over all triplets (i,j,k) with the vertex at site 
j , i and k are nearest neighbors of j . Note that H^i,^ con- 
tains an angular term which is a sort of angular stiffness 
which is essential to stabilizing the diamond lattice (in 
fact, assigning the bonds length alone is not sufficient, 
because the lattice has 3A^ translational degrees of free- 
dom, and the bond length only imposes 27V constraints). 
A single MC step is performed in the following manner: 
an atom of species St at position is randomly chosen 
and a transition to the state r^, 5'^ is attempted. The 
change in energy is then calculated and the move is ac- 
cepted or rejected according to the standard Metropolis 
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criterion. Note that not necessarily both and S'^ have 
to be different from the initial values r^, Si. We sim- 
ulated systems of sizes L = 4 (iV=512), 6 (iV=1728), 8 
(iV=4096), 10 (iV=8000), 12 (iV=13824), 14 (A^^21952), 
and 18 (A^=46656). Periodic boundary conditions were 
used. A number ranging from 5 x 10"* MCS for the small- 
est systems to 5 x 10^ MCS were discarded to thermalize 
the system. The typical number of MC steps for sam- 
pling ranges from 2 x 10^ to 6 x 10^. For each system 
size we performed 10-50 independent runs, using differ- 
ent random number sequences to achieve a satisfactory 
statistical error on the averages of the sampled quanti- 
ties. During a simulation the volume is kept constant at 
a value corresponding to the lattice constant oq of the 
ferromagnet, for consistency to the simulations of Refs. 
127113 In Simulations were performed at fixed chemical po- 
tential /LtB = fJ-o varying the temperature T, as well as at 
fixed T = Tq varying /ib (see Sec. Illlll . In this study we 
deal with antiferromagnetic ordering and consider atoms 
sitting on an "even" site to belong to a sublatticc SLi and 
those sitting on an "odd" site to belong to a different sub- 
lattice SL2. During the simulations we sampled the SW 
energy, and the fraction of A and B particles in SLi and 
SL2. Using these quantities as input, we calculated all 
the thermodynamic quantities needed by employing the 
histogram reweighting method^ 

The order parameter M is defined by the absolute 
value of the staggered magnetization m"*", i.e. 



and 



M = (\m^ 



SL2 



(5) 



(6) 



where the first and second sums in Eq. Q are performed 
over spins belonging to the SLi and SL2 sublattices, re- 
spectively. The finite lattice staggered susceptibility of 
the order parameter is 



X+^N{{m+^)-{\m+\f)/kBT, 



(7) 



where is the Boltzmann constant and T the absolute 
temperature; the staggered susceptibility is 



X+ ^ N {{{m+f) ^ {m+f) /keT. 



(8) 



We also determined the average concentration of the B 
species 



(9) 



where ns is the number of B particles in the system; the 
reduced fourth order cumulant of the order parameter 



C/4-1 



3(m+2)2^ 



(10) 



the specific heat 

a = N{kBTf {{E^) - (Ef) 



(11) 



where E is the total energy. In our analysis we consid- 
ered, in addition, the logarithmic derivatives of and 
of U4 with respect to T. To calculate them we used the 
relation 

Moreover, the derivatives of a thermodynamic quantity 
X with respect to T were calculated using its cross- 
correlation with the energy, i.e. 



d\X\ 
dT 



{\X\E)-{\X\){E). 



(13) 



III. MONTE CARLO RESULTS 

A. Phase diagram 

We started the present study with a rough determina- 
tion of the phase diagram in the (/iB, T) plane. We per- 
formed simulations at different values of T ranging from 
0.05 eV to 0.35 eV.^ At fixed T we swept the chemical 
potential /iB from to 4 eV at intervals of 0.1 eV, while 
/iA was kept constant at 1 eV. For each value of T we de- 
termined the value of /iB at which the maximum value of 
X occurs (see Fig. ^a)). The behavior of concentration 
cb as the chemical potential is swept is shown in Fig. 
^b) . The solid dots indicate the locations of the peaks 
in the finite lattice ordering susceptibility. The result- 
ing phase boundary is plotted in Fig. |2Ia), whereas Fig. 
|2Ib) shows the phase diagram in the (cb,T) plane. Note 
that this procedure gives only an approximate phase di- 
agram, which does not take into account extrapolation 
to the thermodynamic limit. It is, however, very use- 
ful as it provides us with information on where to focus 
in the (/ZB,r) plane with higher resolution simulations 
and with finite size scaling analysis. A transition point 
from the disordered to the ordered phase is estimated at 
/iB ~ 1.42 eV and T ~ 0.34 eV. We decided, therefore, to 
keep /iB fixed at the above value and to run simulations 
in a neighbor of T 0.34 eV. This corresponds to mov- 
ing along the y-axis, i.e. perpendicularly to the phase 
diagram of Fig. \^&). We also performed simulations 
moving parallel to the phase diagram at T = 0.1 eV and 
varying /iB. The results, not shown in this paper, are 
consistent with those obtained by moving in the orthog- 
onal direction. Note the slight asymmetry of the curves 
in Fig. [5]which is a consequence of the three body inter- 
actions. The phase boundary shows no hysteresis, and a 
finite size scaling analysis (descibed in the next section) 
indicates that the transition is 2nd order. 
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FIG. 1: Typical data used to determine the critical points, 
(a) Plot of X/i vs /iB. (b) Concentration of the B species 
vs fiB- The bold circles show the location of the transitions. 
L = 4 and the temperature is T = 0.1 eV in both plots. The 
error bars are smaller than the size of symbols. 



B. Critical behavior 

The critical exponent ly can be determined indepen- 
dently of any other critical quantity, and therefore more 
accurately. It can be shown that the slope of the cumu- 
lant at the critical temperature Tc, or at any other 
point in the critical region, apart from higher order cor- 
rections, scales with system size as L^^'^M' The same 
scaling behavior is exhibited by the logarithmic deriva- 
tive of any power of the staggered magnetization. We 
calculated the derivative of U4, of InM and of ln(m"'"^) 
with respect to T, using Eq. (fT^ . Figure |21 displays a 
log-log plot of the maximum of these derivatives vs. L. 
As expected, they show a linear behavior and they are 
to a very good approximation parallel. No indication of 
correction to scaling are evident. The slopes found af- 
ter a linear fit to the data are 1.579 ± 0.056, 1.619 ± 
0.031, 1.617 ± 0.028, respectively. The weighted average 
of these values gives v = 0.620 ± 0.008, which, within 
the error bars, is in reasonable agreement with the Ising 
value 0.6295 ± 0.0009 found by MC simulations^ but 
is quite different from the theoretically expected Fisher 
renormalized ly' = iy/{l — a) ~ 0.702. 

Various thermodynamic quantities exhibit an extreme 
at a certain temperature Tc{L), which depends strongly 



FIG. 2: Phase diagram of the AF compressible Ising model at 
constant volume in the {ij,b,T) plane (a), and in the (cb,T) 
plane (b). The second order line separates the ordered- 
disordered phases. Estimates are for L = 4. Note the slight 
asymmetry of the curve which reflects the asymmetry of our 
model. The error bars if not visible are smaller than the size 
of symbols. 

on that quantity and on the system size. In the asymp- 
totic regime of large systems the following scaling law 
holds 

Te(L) ^ + A,L-l/^ (14) 

where the subscript x indicates that the prefactor A de- 
pends on the quantity considered. Once v is determined, 
Eq. (|14(l enables us to extrapolate Tc{L) to the ther- 
modynamic limit L ^ 00. In Fig. ^ we have plotted 
Tc{L) of the various quantities examined against L^^l'^ . 
The data follow very close the linear behavior expected, 
except the two smaller lattice sizes L = 4, 6 for which 
corrections to scaling are required. We have therefore 
excluded those data from the extrapolations. The full 
lines in Fig. 0] are linear fits that account for both the 
errors on x and y coordinates. The intercepts on the 
T-axis (L = 00) are very close to each other, however, 
to account for the slight deviations, we have weight av- 
eraged the extrapolated values. The final estimation is 
Tc = 0.34404 ±0.00006 eV. 

The critical exponents 7/1^ and lijv can be directly 
determined from the finite scaling of and M, respec- 
tively at Tc. In the asymptotic regime these quantities 
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FIG. 3: Determination of the critical exponent v from the size 
dependence of the maxima of thermodynamic response func- 
tions. From top to bottom: ain(m+^>/ar (D2), d\uM/dT 
(Dl), dUi/dT (DU4). The error bars are smaller than the 
size of symbols. 
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FIG. 5: Determination of the critical exponents by finite scal- 
ing relations, (a) Linear fit of X^(Tc) (Eq. 1151 '). which pro- 
vides ') /v- (b) Linear fit of M(Tc) (Eq. lltjll '). which provides 
j3/v. Full lines are fits, dotted lines are just extensions of the 
full ones. The error bars if not visible are smaller than the 
size of symbols. 



FIG. 4: Extrapolations of the finite system "critical tempera- 
ture"to the thermodynamic limit for /is = 1.42 eV. From top 
to bottom: dUi/dT, dln{m+^)/dT, dlnPI/dT, x, dM/dT, 
and Cv Full lines are fits of Eqs. 1141 to the corresponding 
data. Dotted lines are just extensions of the full ones. The 
error bars if not visible are smaller than the size of symbols. 

scale as 

X+ ~ L^/^ (15) 
M - L-^^/". (16) 

Figure |5l shows log-log plots of the data. If, again, we ex- 
clude the lowest lattice size, the data are found to follow 
a straight line very well. From a linear fit to the sus- 
ceptibility we get j/iy = 2.017 ± 0.041. Using the value 
previously found for i^, we get 7 = 1.25 ± 0.04, which, 
within the errors, is in agreement with the Ising value 
7ising = 1.2390 ± 0.0071 determined by e-expansion,^ 
and 1.237 ± 0.002 by MC simulations&2^ Analogously, 
we find /3/i/ = 0.491 ± 0.057, or /3 = 0.305 ± 0.039, 
which, within the errors, agrees with the Ising value 
/3i,i„g = 0.3270 ± 0.00154 as well. We have also tried 
to determine the critical exponent of the specific heat a. 



This is however, more difficult to measure because of the 
presence of an additional fitting parameter. The specific 
heat is, indeed, expected to scale at Tc as 

a-Si+SaL"/^ (17) 

The fit of Bi, B2, and a/v, not shown here, gives ajv — 
0.28 ± 0.09, or a = 0.17 ± 0.06, which is in reason- 
able agreement with the Ising value aising = 0.110 ± 

002.35'36 

The fourth order cumulant J74 is an important quan- 
tity to determine the kind of a phase transition and also 
to provide an independent determination of the critical 
temperature. The curves U4,l{T) plotted for different L 
vs T for large L all cross at Tc.^'* Moreover, the value 
C/4 (Tc) strongly depends on the kind of transition. It has 
been found that U4 ~ 0.27052 for the mean-field univer- 
sality class^, U4 c± 0.47 for the three-dimensional Ising 
model, fii and U4 ~ 0.5 for a first order transition. If the 
asymptotic regime has not yet entered, however, curves 
with different L will cross at different points. Neverthe- 
less, it is still possible to extrapolate the crossing point to 
L — > 00. The procedure is described in Ref. 0. Figure El 
displays the behavior of U4 for our system. For L >8 the 
different curves cross with very good approximation at 
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FIG. 6: Plot of the fourth-order cumulant Ui vs L. The 
horizontal dotted line indicates the crossing point value of 
the rigid Ising model. The full lines are just a guide for the 
eyes. The error bars if not visible are smaller than the size of 
symbols. 



the Ising value. The value obtained for the critical tem- 
perature using this method is consistent with the value 
previously determined. 

For a d-dimensional system for which the hyperscaling 
relation dv = 2 — a is valid, the following finite-size laws 
hold in the vicinity of the critical point 



(18) 
(19) 



where t = 1 — T/Tc- In a scaling plot of Xfi^ '^^^ (resp. 

ML'^I'') vs |1 -T/TclL^/'' one should, therefore, observe 
a collapsing of the data. This is exactly what we found, 
as Fig. [3 demonstrates, and is a further evidence of the 
consistency of the critical exponents and Tc previously 
determined. 



interactions^. The reasons of these disagreements are 
not clear and should be further investigated. Needless 
to say, however, that our conclusions should be viewed 
within the context of any numerical work, keeping in 
mind the finiteness of the systems used in the simula- 
tions. It is, therefore, clear that, on the basis of these 
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FIG. 7: Data collapsing of the rescaled susceptibility X/; (a), 
and of the order parameter M (b) vs rescaled temperature. 
The error bars if not visible are smaller than the size of sym- 
bols. 



IV. CONCLUSIONS 

We have performed Monte Carlo simulations of the 
compressible Ising model with antiferromagnetic inter- 
actions under constant volume conditions, in the semi- 
grand-canonical ensemble. Elastic forces are included by 
the Stillinger- Weber potential. The behavior of all crit- 
ical quantities analyzed strongly indicated the presence 
of a closed second order line with the critical exponents 
of the (rigid) Ising model. This is in contrast with theo- 
ries as they predict the occurrence of Fisher renormalized 
exponents. Disagreement was also found in the simula- 
tions of exactly the same model but with ferromagnetic 



data, the occurrence of a (slow) crossover toward Fisher 
renormalized exponents cannot be completely ruled out. 
A deeper investigation of this issue would require simu- 
lations on much larger system sizes, which is, basically, 
unfeasible with the present computer power. It would 
be, however, very interesting to investigate the critical 
behavior under stronger coupling conditions. 
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